High Temporal Resolution Estimation of Hemodynamic Response from Event-related fMRI

Introduction Estimation of hemodynamic response associated with brain activation following various types of events, i.e., motor, sensory, cognitive, often relies on the solving of a deconvolution problem. Deconvolution methods employed to date have been limited in their ability to investigate details of the hemodynamic response function (HRF) due to the limited temporal resolution of the deconvolved HRF [e.g., 1-4]. In particular, these methods provided estimates with a temporal resolution that is fixed at the volume acquisition time, TR. In addition, they could not provide estimates that were truly slice-specific and needed to invoke an interpolation to correct for what is known as the slice timing problem [5]. We present here an approach to computation that automatically ensures a slice-specific solution and allows for a solution at a temporal resolution greater than one TR. This method is demonstrated with an experimental event-related fMRI data set. Methods We assume a general linear model approach. The convolution can be expressed as f(t) = g(t)*h(t) + B(t) + η(t), where f(t) is the response that is measured, h(t) is the HRF, g(t) is a stimulus or event function comprised of delta functions, B(t) is a function describing the fMRI baseline, and η(t) is attendant noise. Presentation of the stimuli should be jittered in time [6,7] to aid in optimal estimation. Explicit writing out of the equation for f(t) in the discrete case yields f = Xh + n, where f is a vector of the fMRI time series, X is the experimental design matrix, and n is the noise vector. The rows of the design matrix consist of the coefficients from the simultaneous linear equations that result from the convolution, plus terms used to fit B(t). The general linear model provides a least squares estimate [6] which, assuming uncorrelated normally distributed noise, can be written as h= (XX)Xf. The rows for the X matrix arise from the set of linear simultaneous equations that result from writing out the convolution in discrete form: f(t) = Σ g(τ) h(t-τ)δ. Many deconvolution schemes adopt a unit of time, δ, for discretizing f, g, and h that is one TR, since fMRI time series are measured one datum per TR. Our method adopts a smaller time unit that is equal to TR divided by the number of axial slices acquired; this approach tacitly assumes that the EPI acquisition style acquires slices that are evenly spaced throughout the TR. Use of this smaller division of time ensures that each row of the X matrix describes what is occurring for one slice at a time. From the X matrix, one can then construct a smaller submatrix that contains only those rows that corresponds to a specific slice. This smaller slice-specific matrix can then be used to solve for h, as indicated above. Since δ can be on the order of 100ms, the number of the δ units needed to span the interval of a HRF could total over 300. To reduce this number to a computationally manageable size we adopt the following approximation. We select a HRF interval that is spanned by a number of δs, q, that is divisible by a small power of two or a product of two small prime numbers, .e.g., 8. The interval is then partitioned into q/8 sections. The values of h in each of these sections are then approximated as proportions of the ‘distance’ between the endpoints of the sections; thus each section is approximated by a straight line. The values of h between the endpoints are then substituted back into the original simultaneous equations. Like terms for the remaining h values can be grouped and the X matrix is then constructed with fewer unknowns. The number that is chosen to divide q can be picked as small as the signal-to-noise ratio of the experiment will allow; the largest this number would logically be would the number of slices. From the X matrix one can compute a t-statistic [8] for each value of h. We obtained data on a healthy volunteer using our standard fMRI experimental parameters: Echo planar images were recorded at 3T (Siemens Tim Trio) with TE/TR/flip=29ms/2800ms/80, 31 interleaved axial slices, matrix=128x128, and 256mm x 256mm FOV. The volunteer viewed in the scanner a display that was fundamentally a black background with a small white plus sign in the center. At jittered times a flashing checkerboard (black and white, 10Hz), with the same plus sign centered in it, was presented for 1.0 second. The inter-stimulus interval ranged from 5 to 15 seconds. To minimize head motion a bitebar was used. Before deconvolution, the data set underwent motion correction, and a Hamming spatial filter was applied to improve SNR [9]. To minimize the effects of non-equilibrium spin history the first four volumes of the data set were removed. No other treatment of the data was performed. Results and Discussion Following deconvolution of the experimental data set with our method, large regions of activation were found around the calcarine fissure. Figure 1 shows the activation appearing in two adjoining slices. Regions of interest were drawn for the two slices around the activated voxels. A plot for the ROI-averaged HRFs for the slices is given in Figure 2. These curves were calculated at a temporal resolution of 1445ms; the value of q was 320 and the power-of-two number was 16. The number of voxels with t > 3.5 is given. Because of the interleaved slice acquisition of the experiment, the adjoining slices were acquired ~1.4 seconds apart. Thus, it is remarkable how well the two curves agree with each other in location of maxima, main peak width and shape. The SNR was high enough to allow deconvolution at higher temporal resolution. Figure 3 shows the HRF resulting from averaging over 362 voxels in four contiguous slices; the power-of-two number was 8 and the resulting temporal resolution is ~723ms. We emphasize that Figs. 2 and 3 were extracted from the same data resulting from one scan measured with a TR of 2.8s.